<!-- XHTML 1.0 Strict -->
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<meta name="author" content="Rachid Touzani" />
<meta name="keywords" content=" css, dropdowns, dropdown menu, drop-down, menu, navigation, nav, horizontal, vertical left-to-right, vertical right-to-left, horizontal linear, horizontal upwards, cross browser, internet explorer, ie, firefox, safari, opera, browser, lwis" />
<meta name="description" content="Clean, standards-friendly, modular framework for dropdown menus" />
<link href="../css/dropdown.vertical.css" media="screen" rel="stylesheet" type="text/css" />
<link href="../css/default.ultimate.css" media="screen" rel="stylesheet" type="text/css" />

<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<title>Fluid Dynamics Demo 1</title>
<link rel="stylesheet" type="text/css" href="../doxygen.css" />
<link rel="stylesheet" type="text/css" href="../tabs.css"  />
</head>

<body bgcolor="#FFFFFF" link="#FF0000" vlink="#FF00FF" alink="#FF00FF">
<p align="center"><a href="../index.html"><img src="../im/ofeli.gif" width="150" border="0"></a></p>

<center>
  <div id="cse" style="width: 60%;">Loading</div>
  <script src="http://www.google.com/jsapi" type="text/javascript"></script>
  <script type="text/javascript"> 
     google.load('search', '1', {language : 'en', style : google.loader.themes.SHINY});
     google.setOnLoadCallback(function() {
        var customSearchOptions = {};  var customSearchControl = new google.search.CustomSearchControl(
        '012396140824982761142:-qrywxcfi_o', customSearchOptions);
        customSearchControl.setResultSetSize(google.search.Search.FILTERED_CSE_RESULTSET);
        customSearchControl.draw('cse');
     }, true);
   </script>
   <style type="text/css">
     .gsc-control-cse {
     font-family: Verdana, sans-serif;
     border-color: #DAE0E5;
     background-color: #DAE0E5;
   }
   .gsc-control-cse .gsc-table-result { font-family: Verdana, sans-serif; }
   input.gsc-input { border-color: #B6BEC5; }
   input.gsc-search-button {
     border-color: #B6BEC5;
     background-color: #D0D1D4;
   }
   .gsc-tabHeader.gsc-tabhInactive {
     border-color: #999999;
     background-color: #EEEEEE;
   }
   .gsc-tabHeader.gsc-tabhActive {
     border-color: #999999;
     background-color: #999999;
   }
   .gsc-tabsArea { border-color: #999999; }
   .gsc-webResult.gsc-result,
   .gsc-results .gsc-imageResult {
     border-color: #FFFFFF;
     background-color: #FFFFFF;
   }
   .gsc-webResult.gsc-result:hover,
   .gsc-imageResult:hover {
     border-color: #D2D6DC;
     background-color: #EDEDED;
   }
   .gsc-webResult.gsc-result.gsc-promotion:hover {
     border-color: #D2D6DC;
     background-color: #EDEDED;
   }
   .gs-webResult.gs-result a.gs-title:link,
   .gs-webResult.gs-result a.gs-title:link b,
   .gs-imageResult a.gs-title:link,
   .gs-imageResult a.gs-title:link b { color: #0568CD; }
   .gs-webResult.gs-result a.gs-title:visited,
   .gs-webResult.gs-result a.gs-title:visited b,
   .gs-imageResult a.gs-title:visited,
   .gs-imageResult a.gs-title:visited b { color: #0568CD; }
   .gs-webResult.gs-result a.gs-title:hover,
   .gs-webResult.gs-result a.gs-title:hover b,
   .gs-imageResult a.gs-title:hover,
   .gs-imageResult a.gs-title:hover b { color: #0568CD; }
   .gs-webResult.gs-result a.gs-title:active,
   .gs-webResult.gs-result a.gs-title:active b,
   .gs-imageResult a.gs-title:active,
   .gs-imageResult a.gs-title:active b { color: #0568CD; }
   .gsc-cursor-page { color: #0568CD; }
   a.gsc-trailing-more-results:link { color: #0568CD; }
   .gs-webResult .gs-snippet,
   .gs-imageResult .gs-snippet,
   .gs-fileFormatType { color: #5F6A73; }
   .gs-webResult div.gs-visibleUrl,
   .gs-imageResult div.gs-visibleUrl { color: #5F6A73; }
   .gs-webResult div.gs-visibleUrl-short { color: #5F6A73; }
   .gs-webResult div.gs-visibleUrl-short { display: none; }
   .gs-webResult div.gs-visibleUrl-long { display: block; }
   .gs-promotion div.gs-visibleUrl-short { display: none; }
   .gs-promotion div.gs-visibleUrl-long { display: block; }
   .gsc-cursor-box { border-color: #FFFFFF; }
   .gsc-results .gsc-cursor-box .gsc-cursor-page {
      border-color: #999999;
      background-color: #FFFFFF;
      color: #0568CD;
   }
   .gsc-results .gsc-cursor-box .gsc-cursor-current-page {
      border-color: #999999;
      background-color: #999999;
      color: #0568CD;
   }
   .gsc-webResult.gsc-result.gsc-promotion {
      border-color: #D2D6DC;
      background-color: #D0D1D4;
   }
   .gsc-completion-title { color: #0568CD; }
   .gsc-completion-snippet { color: #5F6A73; }
   .gs-promotion a.gs-title:link,
   .gs-promotion a.gs-title:link *,
   .gs-promotion .gs-snippet a:link { color: #0066CC; }
   .gs-promotion a.gs-title:visited,
   .gs-promotion a.gs-title:visited *,
   .gs-promotion .gs-snippet a:visited { color: #0066CC; }
   .gs-promotion a.gs-title:hover,
   .gs-promotion a.gs-title:hover *,
   .gs-promotion .gs-snippet a:hover { color: #0066CC; }
   .gs-promotion a.gs-title:active,
   .gs-promotion a.gs-title:active *,
   .gs-promotion .gs-snippet a:active { color: #0066CC; }
   .gs-promotion .gs-snippet,
   .gs-promotion .gs-title .gs-promotion-title-right,
   .gs-promotion .gs-title .gs-promotion-title-right *  { color: #333333; }
   .gs-promotion .gs-visibleUrl,
   .gs-promotion .gs-visibleUrl-short { color: #5F6A73; }
   </style>
</center>
</p>
&nbsp;

<div id="mainmenu">
<div class="text">

 <div class="tabs">
   <ul class="tablist">
   <li><a href="../index.html"><span>Home</span></a>
   <li><a href="../overview.html"><span>Overview</span></a>
   <li><a href="../html/index.html"><span>Class Documentation</span></a>
   <li><a href="../fformats.html"><span>File Formats</span></a>
   <li><a href="../tutorial.html"><span>Tutorial</span></a>
   <li class="current"><a href="../demos.html"><span>Demo Codes</span></a>
 </div>

 </div>
 </div>

&nbsp;
<p>
<h2>Fluid Dynamics Demo 1: A 2-D incompressible CFD problem using quadrilateral elements</h2>

<p><table border="2" cellpadding="2" cellspacing="0" >
<tr>
<td align="center" width="35"><a href="truss.html"><img src="../im/backward.gif" border="0"></a></td>
<td align="center" width="35"><a href="../demos.html"><img src="../im/top.gif" border="0"></a></td>
<td align="center" width="35"><a href="tins2.html"><img src="../im/forward.gif" border="0"></a></td>
</tr>
</table>


&nbsp;&nbsp;
<h1>The main program</h1>

<p>
<SPAN class=TEXT>
Let us present an <span class=LOGO>OFELI</span> program that solves the 2-D time dependent incompressible 
Navier-Stokes equations.
</SPAN>

<ul type=square>
<li><SPAN class=TEXT>We start by including the classes for fluid dynamics. We also
use here a user defined class for prescibing data:
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>   #include "OFELI.h"
   #include "Fluid.h"
   #include "User.h"
</font></pre>
</td></tr></table>
</SPAN></p>

<li><SPAN class=TEXT>We define an instance of class <span class=VAR>
<a href="../html/classOFELI_1_1IPF.html">IPF</a></span> to deal with project data.
We extract useful information from this instance. Note that all these parameters were
described in the previous examples.
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>   IPF data("tiff2 - 1.0",argv[1]);
   int output_flag = data.getOutput();
   int save_flag = data.getSave();
   int pres_flag = data.getIntPar(1);
   theTimeStep = data.getTimeStep();
   theFinalTime = data.getMaxTime();
</font></pre>
</td></tr></table>
</SPAN></p>

<li><SPAN class=TEXT>As usual we construct a 
<span class=VAR><a href="../html/classOFELI_1_1Mesh.html">Mesh</a></span> instance.
In addition, we present here another method for prescribing problem data. This one
consists in constructing a class that inherits from abstract class
<span class=var><a href="../html/classOFELI_1_1UserData.html">UserData</a></span>.
This class is presented after the main program.
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>   Mesh ms(data.getMeshFile());
   User ud(ms);
</font></pre>
</td></tr></table>
</SPAN></p>

<li><SPAN class=TEXT>
The linear system uses a symmetric matrix. For this, we
we use an instance of class <span class=VAR><a href="../html/classOFELI_1_1SkSMatrix.html">SkSMatrix&lt;double&gt;</a></span>.
Then we declare a vector <span class=VAR>b</span> that will contain the right-hand side
and then the solution. We also need a vector <span class=VAR>u</span> that contains
the solution at the previous time step:
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>   SkSMatrix&lt;double&gt; A(ms);
   Vect&lt;double&gt; b(ms), u(ms);
</font></pre>
</td></tr></table>
</SPAN></p>

<li><SPAN class=TEXT>
We now use the instance <span class=VAR>ud</span> to prescribe initial solution,
Dirichlet boundary condition, body force and source force (tractions):
<span class=VAR><a href="../html/classOFELI_1_1SkSMatrix.html">SkSMatrix&lt;double&gt;</a></span>.
Then we declare a vector that will contain the right-hand side and then
the solution:
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>   ud.setInitialData(u);
   Vect&lt;double&gt; bc(ms);
   ud.setDBC(bc);
   Vect&lt;double&gt; body_f(ms);
   ud.setBodyForce(body_f);
   Vect&lt;double&gt; bound_f(ms);
   ud.setSurfaceForce(bound_f);
</font></pre>
</td></tr></table>
</SPAN></p>

<li><SPAN class=TEXT>In order to store the solution we define 2 instances of class
<span class=VAR><a href="../html/classOFELI_1_1IOField.html">IOField</a></span>
to store respectively the velocity and the pressure. The names of the storage
files were primarily obtained from the project file.
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>   IOField vf(data.getMeshFile(),data.getString("v_file"),ms,IOField::OUT);
   IOField pf(data.getMeshFile(),data.getString("p_file"),ms,IOField::OUT);
</font></pre>
</td></tr></table>
</SPAN>
</font></pre>
</td></tr></table>
</SPAN></p>

<li><SPAN class=TEXT>
Once all the data are read, we start the time loop. At first, we set to <span class=TT>0</span>
the right-hand side that will be assembled.
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>   TimeLoop {
      b = 0;
</font></pre>
</td></tr></table>
</SPAN>
</font></pre>
</td></tr></table>
</SPAN></p>

<li><SPAN class=TEXT>The global matrix and right-hand side are, as usual, filled
by looping over elements. For each element, an instance of class 
<span class=VAR><a href="../html/classOFELI_1_1NSP2DQ41.html">NSP2DQ41</a></span>
is considered. This class deals with quadrilateral elements. In this formulation, a penalty
term takes into account the divergence free constraint. We add for this the function 
<span class=VAR>Penal</span> that adds the contribution of penalty term with a penalty
parameter equal here to <span class=TT>10<sup>7</sup></span>. We also add the (lumped) mass
 matrix and the viscosity term with viscosity equal here to <span class=TT>0.1</span>.
We then add the body forces using the <span class=VAR>User</span> instance.
Assembly of element matrix is performed only at the first time step, the matrix being
constant. Assembly of the element right-hand side is also made.
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>      MeshElements(ms) {
         NSP2DQ41 eq(theElement,u,theTime);
         eq.LMass(1./theTimeStep);
         eq.Penal(1.e07);
         eq.Viscous(0.1);
         eq.RHS_Convection();
         eq.BodyRHS(ud);
         if (theStep==1)
            eq.ElementAssembly(A);
         eq.ElementAssembly(b);
      }
</font></pre>
</td></tr></table>
</SPAN></p>

<li><SPAN class=TEXT>Since we want to take account for boundary forces, we have
to loop over sides (here edges) and construct element loads. This has also been 
considered in previous lessons.
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>      MeshSides(ms) {
         NSP2DQ41 eq(theSide);
         eq.BoundaryRHS(ud);
         eq.SideAssembly(b);
      }
</font></pre>
</td></tr></table>
</SPAN></p>

<li><SPAN class=TEXT>We can now prescribe boundary conditions by penalty method
and solve the resulting linear system by factorization and backsubstitution.
We finally assign the resulting solution to the vector <span class=VAR>u</span>.

<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>      A.Prescribe(ms,b,bc,theStep-1);
      A.solve(b);
      u = b;
</font></pre>
</td></tr></table>
</SPAN></p>

<li type=square><SPAN class=TEXT>We can save the resulting velocity vector
for plotting or any other purpose. For this, the class
<span class=VAR><a href="../html/classOFELI_1_1IOField.html">IOField</a></span>
is used.<br> Note that we need to convert to class 
<span class=VAR><a href="../html/classOFELI_1_1Vect.html">Vect&lt;double&gt;</a></span>
in order to save the vector with some mesh information.
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>      u.setTime(theTime);
      if (output_flag > 0)
         cout << u;
      if (save_flag)
         v_file.put(u);
</font></pre>
</td></tr></table>
</SPAN></p>

<li><SPAN class=TEXT>If the flag <span class=VAR>pres_flag</span> is true, then
we will compute <i>a posteriori</i> the pressure field. This is easy when
the penalty formulation is used: We start by computing a piecewise constant pressure field, then use a smoothing procedure to retrieve a continuous piecewise 
linear field.
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>      if (pres_flag) {
         double pres = 0.;
         Vect<double> pm(ms,1,NODE_DOF), ep(ms,1,ELEMENT_DOF), p(ms,1,NODE_DOF);
         p.setTime(theTime);
         p.setName("Pressure");
         Reconstruction rr(ms);
         MeshElements(ms) {
            NSP2DQ41 eq(theElement,u,theTime);
            pres += eq.Pressure(1.e07);
         }
         ep(theElementLabel) += pres;
         rr.P0toP1(ep,p);
         if (output_flag > 0)
            cout << p;
         if (save_flag)
            p_file.put(p);
      }
</font></pre>
</td></tr></table>
</SPAN></p>
</ul>

&nbsp;&nbsp;
<h1>The class 'User'</h1>

<SPAN class=TEXT>
We define here the class <span class=VAR>User</span> that defined the Dirichlet boundary
conditions. All other data (body, surface forces and initial conditions) are assumed
null, so they don't need be defined. This class inherits from the abstract class
<span class=VAR><a href="../html/classOFELI_1_1UserData.html">UserData</a></span>.
The constructor using mesh is already implemented in
<span class=VAR><a href="../html/classOFELI_1_1UserData.html">UserData</a></span>.
We need here to define the function <span class=VAR>BoundaryCondition</span> 
in a very simple way.

<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>class User : public UserData<double>
{
 public:

   User(Mesh &mesh) : UserData<double>(mesh) { }

   double BoundaryCondition(const Point<double> &x, int code, double time=0., size_t dof=1)
   {
      time = 0;
      dof = 1;
      double ret = 0.0;
      if (code==1)
         ret = 0.0;
      if (code==2)
         ret = 1.;
      return ret;
   }
};
</font></pre>
</td></tr></table></p>

&nbsp;&nbsp;
<h1>An example</h1>

<SPAN class=TEXT>
Let us run this program with the data presented in the following project file:
<p><table bgcolor="#e0f8f7" border="1" frame=box rules=none width=700>
<tr><td>
<pre><font color="blue" size=3>&lt;?xml version="1.0" encoding="ISO-8859-1" ?&gt;
&lt;OFELI_File&gt;
&lt;info&gt;
   &lt;title&gt;&lt;/title&gt;
   &lt;date&gt;&lt;/date&gt;
   &lt;author&gt;&lt;/author&gt;
&lt;/info&gt;
&lt;Project name="cavity"&gt;
   &lt;time_step value="0.1"/&gt;
   &lt;max_time value="1.0"/&gt;
   &lt;verbose value="1"/&gt;
   &lt;output value="1"/&gt;
   &lt;mesh_file value="cavity.m"/&gt;
   &lt;parameter label="v_file" value="cavity.v"/&gt;
   &lt;parameter label="p_file" value="cavity.p"/&gt;
&lt;/Project&gt;
&lt;/OFELI_File&gt;</font></pre>
</td></tr></table></p>

<p>
Note that only execution and time integration parameters are given.
All other data were given through the class <span class=var>User</span>.</p>

<p><table align="center" border="2" cellpadding="2" cellspacing="0" >
<tr>
<td align="center" width="35"><a href="truss.html"><img src="../im/backward.gif" border="0"></a></td>
<td align="center" width="35"><a href="../demos.html"><img src="../im/top.gif" border="0"></a></td>
<td align="center" width="35"><a href="tins2.html"><img src="../im/forward.gif" border="0"></a></td>
</tr>
</table>

&nbsp;
<div id="foot_bar">
        Copyright &copy; 1998-2018 Rachid Touzani&nbsp;
</div>

</body>
</html>
